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Hamiltonian optics of hyperbolic polaritons in nanogranules 
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Semiclassical quantization rules and numerical calculations are applied to study polariton modes 
of materials whose permittivity tensor has principal values of opposite sign (so-called hyperbolic 
materials). The spectra of volume- and surface-confined polaritons are computed for spheroidal 
nanogranules of hexagonal boron nitride, a natural hyperbolic crystal. The field distribution created 
by polaritons excited by an external dipole source is predicted to exhibit ray-like patterns due to 
classical periodic orbits. Near-field infrared imaging and Purcell-factor measurements are suggested 
to test these predictions. 


Introduction. Recently much interest has been at¬ 
tracted to a class of uniaxial materials whose axial ^ and 
tangential £± permittivities have opposite signs. These 
hyperbolic materials (HM) possess extraordinary rays 
with unusual properties. In this Letter we focus on po¬ 
lar dielectric HM 1 5 where the extraordinary rays are 
phonon-polariton collective modes. Our results may also 
apply to other HM, including ferromagnets 6 , magnetized 
plasmas 7 , artificial metamaterials 8 , layered superconduc¬ 
tors 9,10 , and liquid crystals 11 . 

The basic properties of hyperbolic polaritons are as 
follows. Their isofrequency surfaces o;(p) = const in mo¬ 
mentum space p = (Px,Py,Pz) are hyperboloids. In a 
broad range of |p| from the free-space photon momentum 
uj/c to an upper cutoff imposed by microscopic structure, 
these hyperboloids can be approximated by cones [Fig¬ 
ure 1(a)] 

Hb{ P,w) = e z (uj)p 2 z +£±{lu)(p 2 x +p 2 y ) = 0. (1) 

The group velocity v(p) = d p u: is always orthogonal 
to the isofrequency surface. Hence, within the conical 

approximation it has a fixed angle a = tan -1 

with respect to the optical axis. Such a strictly direc¬ 
tional propagation of polaritons may be used for sub- 
diffractional focusing 12,13 and super-resolution imaging 
known as ‘hyperlensing’ 8,14-16 . Since the high-momenta 
polaritons remain immune to evanescent decay, volume- 
confinement of polaritons inside nanogranules 3,17,18 is 
possible. Several experimental observations of such 
modes in hexagonal boron nitride (hBN), a natural mid- 
infrared HM, have been reported 2,3,12,13,17,19 . (This lay¬ 
ered insulator is also known to be a premier substrate 20 
or a spacer for van der Waals heterostructures 21,22 .) In 
the far-field spectroscopy 3 the polariton modes of hBN 
nanogranules show up as discrete resonances. Remark¬ 
ably, the spectrum of such resonances was found to de¬ 
pend primarily on the aspect ratio of the granules rather 
than their size or precise shape. Exact solutions 1,6,10 
for spheroidal or spherical shapes enable one to compute 
such spectra but they do not elucidate the underlying 
physical picture. 

In this Letter, we further develop an alternative ray 
optics method 23 that makes connection to the Einstein- 



FIG. 1. (a) A schematic of a polariton isofrequency surface 
and the group velocity v in a HM. The example shown is for 
the case ej_ < 0, > 0. (b) The geometry of the model 

studied. Vector d symbolizes an external dipole source. 


Brillouin-Keller (EBK) quantization 24,25 of a classical 
particle inside a cavity having the same shape as the 
granule. The indefinite permittivity tensor of the HM 
maps on the indefinite Hamiltonian Hb{ p,cj) of the par¬ 
ticle, eq (S5). The EBK quantization rules are valid pro¬ 
vided the classical motion is regular 25 . However, they 
give valuable physical insights even when it is weakly 
chaotic 26,27 or pseudointegrable 28 . Classical motion in 
a spheroidal cavity is completely integrable, and so the 
EBK rules apply directly to spheroidal granules. 

Nanoscale spatial distribution of the electric field pro¬ 
duced by the polariton modes can be measured by scan¬ 
ning near-field optical microscopy 2,13,17 . To model such 
experiments we compute the response of a spheroidal 
nanogranule to an oscillating electric dipole. We calcu¬ 
late the field reflected back to the dipole and the field dis¬ 
tribution it induces in the interior and on the surface of 
the spheroid. Both of them exhibit striking geometrical 
patterns that correspond to periodic orbits of polaritonic 
rays. 

The exact eigenmodes. Consider a granule that has 
a shape of a spheroid with the symmetry axis parallel to 
the optical axis (z- axis) of the permittivity tensor (Fig¬ 
ure lb). We assume the spheroid is prolate, a z > a±. 
(Oblate spheroids can be treated in a similar manner.) 
The cross-section of the granule in the cylindrical coor¬ 
dinates p = a/x 2 + y 2 and z is shown in Figure 2a, c. 
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We define two other sets of coordinates. Outside the 
spheroid, we use the usual spheroidal ones: 

p = asinh^sin# , z = a cosh p cos 6 , ( 2 ) 

where p > p, 0 < 6 < 7 r, a 2 = a 2 — a^_, and 

tanh fj = — . (3) 

a z 

These coordinates are orthogonal and real. Inside, Figure 
2 b and d, we use 

p = — iby/e± sin £ sin 0 , z = cos £ cos 6 , (4) 

where 0 < £ < £, £ < 0 < 7 r — £, and 

tan£ = i—(5) 
a z y/el 

Parameter b = (e^a 2 — 5 j 1 a ^) 1 / 2 is real if £_l < 0, 
> 0 , and imaginary if both the signs are reversed. 
The coordinates (£, 0) are real and leave the permittivity 
tensor diagonal. However, they are nonorthogonal. 

We assume that granule is suspended in vacuum and 
that its size is much smaller than c/uo. In this case the 
quasi-static approximation for the electric field is valid, 
E = — <9 r 4>. The scalar potential 4> is represented by two 
different functions inside and outside the particle: 4>i and 
$ 2 . The inner potential obeys the Walker equation 6 

[s z d 2 z + sj_(dl + d 2 y )\<&i = 0 . ( 6 ) 

The outer potential <£>2 satisfies the Laplace equation. 
The potential and the normal component of the displace¬ 
ment must be continuous across the spheroid surface. In 
the chosen coordinates, these boundary conditions be¬ 
come separable, which enables one to find the analytical 
solutions 1,6 of this eigenproblem: 

oc Pr(cosOPr(cos0)e^ , 0 < £ < £, (7a) 

$2 oc Q™ (cosh 77 ) Py 1 (cos 6)e irn( t > , p > fj . (7b) 

Here m, l are integers, P] n (z) 1 Q] n {z ) are the associated 
Legendre functions of the first and the second kinds, re¬ 
spectively, and (j) is the polar coordinate in the x-y plane. 
The boundary conditions are satisfied provided 1,6 

4 lnP™(cosf) = InQ™(cosh fj) . ( 8 ) 

at, ap 

For each m and l this equation gives us several solutions 
for the eigenfrequency uj (contained implicitly in £ 1 , e z ) 
that can be indexed by another integer n. The total 
number of such solutions is equal to / for m = 0 and 
/ — \m\ + 1 for m ^ 0 (Supporting information). Note 
that eq (SI) depends only on the aspect ratio and not 
the size of the spheroid. This is consistent with the scale- 
invariance of eq ( 6 ). However, the physical picture is 
not clear from this exact solution. Next, we present an 





FIG. 2. The correspondence of the (p,z) and (£, 6) coordi¬ 
nate systems for the spheroid. 1, 2, 3 label the three bound¬ 
ary regions. The green dots separating these regions mark 
the points where the spheroid surface is tangent to the po- 
lariton group velocity v. (a,b) the classically accessible bulk 
region for a wave with a caustic £ = £ c (dotted line) inside 
the spheroid. Lines inside the region are trajectories of wave 
packets, which are straight lines in real space although they 
appear as curves in our choice of coordinates, (c, d) The 
classically accessible boundary region (thick dark line) for a 
surface wave. The caustic (dotted line) extends outside the 
spheroid. 

alternative derivation in terms of a more intuitive ray- 
optics approach. 

Hamiltonian optics. Ray or geometrical optics is a 
well established approach to study propagation of light 
on scales longer than the photon wavelength. HMs are a 
new arena for ray optics in which photons are replaced 
by excitations of much shorter wavelength — polaritons 
- in the case of hyperbolic polar insulator. This ap¬ 
proach has been previously applied to HM of cylindrical 
geometry 23,27 . Here we study a spheroidal granule and 
address both the ray and the wave optics effects within 
the quasi-static approximation. The derivation of the 
ray picture starts with seeking the scalar potential in the 
form 

$ 1 (r) = ^A i (r)e^'W, (9) 

j 

where the phases (or eikonals) Sj (r) vary much faster 
than the amplitudes Aj (r) . Substituting eq (9) into 
eq (6) and keeping only the leading terms, quadratic 
in momenta Pj (r) = d r Sj, one obtains (for each j ) the 
Hamilton-Jacobi equation (S5) of a fictitious classical sys- 
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tem with the ‘optical’ Hamiltonian Hsipj, uo). The EBK 
quantization is possible if the number of different j in 
eq (9) is finite. For the spheroid four terms suffice, cor¬ 
responding to the different sign choices of the momenta. 
To describe 0- and ^-motions one needs two terms each 
because our fictitious particle can propagate in two op¬ 
posite directions between the surface and the caustics. 
The 0-motion has no caustic and one term is enough. 
We must clarify that ‘motion’ and ‘propagation’ refer to 
the geometry of the phase-space flow, not to the actual 
time evolution of coordinates and momenta. The veloc¬ 
ity = d p H b of the fictitious particle deduced from the 
optical Hamiltonian is different from the group velocity 
of an actual polariton 


Olj 


dp 

V / h b = o 


( 10 ) 


However, and v are always parallel to each other. 
Therefore, if a fictitious classical particle with the con¬ 
served energy Hb( p,u;) = 0 and a polariton wavepacket 
of frequency uj are launched at the same initial point 
(p, r), the geometrical shape of their phase-space trajec¬ 
tories will be identical. This identity is well known in the 
Hamiltonian formulation of geometrical optics 29 . Here 
we adopt it for hyperbolic polaritons. The EBK quan¬ 
tization rules 24,30 can also be directly adopted for our 
problem because they are formulated in terms of contour 
integrals in the phase-space. Therefore, to compute po¬ 
lariton eigenmodes of an arbitrary nanostructure made of 
HM, we need to quantize the motion of a single particle 
bouncing inside a cavity of the same shape. 

Two unusual circumstances still have to be handled. 
First, the Hamiltonian of our fictitious particle is indef¬ 
inite. Second, the reflection rule and the corresponding 
phase shift at the surface are determined by the bound¬ 
ary conditions. For spheroidal nanogranule (Figure lb), 
both of these peculiarities prove to be tractable in the co¬ 
ordinate system defined above. Performing the canonical 
transformation to the new momenta p^pe, we find 

rr . Pj-PB _ Vj, 

B sin 2 £ — sin 2 0 sin 2 £ sin 2 9 


The classical motion governed by Hamiltonian eq (11) is 
separable and so integrable. (Unfortunately, in the class 
of smooth convex shapes, only ellipsoids and spheroids 
as their particular case appear to be integrable 31 .) 

The EBK quantization rules are in the form 



7r 
2 


+ 5 = 2,7ns , 


pedO + 25 = — 2tt\ , 


27 TPcf) = 27 Tfl . 


(12a) 

(12b) 

(12c) 


Here each integral is taken over a closed-loop contour in 
a respective coordinate (cf. Figure 2b). The phase shift 


S is constant everywhere on the surface (Supporting in¬ 
formation), as demanded by the separable form of the 
exact solution, eq (7). These integrals can be evaluated 
in terms of elementary functions. Comparing those ex¬ 
pressions with the asymptotic formulas 24 for Legendre 
functions, it is easy to establish the correspondence be¬ 
tween the EBK quantum numbers A, z/, p and the indices 
/, m, n in the exact solution (Supporting information): 

1 = 2z/ + A + /r, ra = /i, n = v . (13) 


We refer to the v > 0 eigenmodes as the bulk modes and 
to those with v = 0 as the surface ones. The scalar po¬ 
tential <f>i of the bulk modes oscillates inside the granule 
along the ‘radial’ direction £ whereas that of the sur¬ 
face modes monotonically increases with £ and reaches a 
maximum at the surface. 

To compare the EBK results with the exact solution, 
we calculated the eigenmode spectra of an hBN spheroid 
as a function of its aspect ratio. The measured 3 optical 
constants of hBN were used except the damping was ne¬ 
glected in order to obtain real solutions for cj. Examples 
of these calculations are shown in Figure 3. The EBK 
is expected to be asymptotically exact at large quantum 
numbers but as one can see from Figure 3, an excellent 
agreement is reached for the bulk modes (the top three 
curves) already for modest /, m, and n. On the other 
hand, the (9, 2, 0) surface mode (the bottom curve) shows 
some deviations from the exact result at intermediate as¬ 
pect ratios. We discuss such modes in more detail below. 

Surface modes. The hyperbolic surface modes 
(HSM) 32 are similar 33 to Dyakonov surface waves 34-36 
of uniaxial materials with positive-definite permittivity 
tensor. However, the HSMs have several new properties. 
Unlike the standard Dyakonov waves, the momenta and 
therefore achievable degree of confinement for the HSM 
are limited only by microscopic (for hBN, atomic) struc¬ 
ture. The HSM are robust to surface defects in the sense 
that there can only be three other fixed directions for the 
defect-scattered wave. This is a stronger angular restric¬ 
tion than the absence of electron backscattering in topo¬ 
logical insulators 37 and graphene 38 . Finally, compared 
to surface plasmons in metals, which lack any direction¬ 
ality, the HSM of polar insulators should exhibit a much 
lower damping as the they are free of electronic losses. 

In the present case of the spheroid, the HSM corre¬ 
spond to the EBK quantum numbers v = 0 (and so to 
n = 0). The ^-coordinate of the caustic is given by 


sin £ c 


m 



A + p + | 


(14) 


see Supporting information. For l and m fixed, £ c is inde¬ 
pendent of the aspect ratio A = a±/a z . The coordinate £ 
of the spheroid surface [eq (3)] increases with A. At large 
enough A, we have £ > £ c , see Figure 2b. This is simi¬ 
lar to bulk modes (i.e., n > 0 modes) except the caustic 
is now very close to the surface. At small A, we have 
£ < £ c , so the caustic extends beyond the surface, Fig¬ 
ure 2d. Momentum p £ is imaginary inside the spheroid 
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Aspect Ratio 

FIG. 3. Eigenfrequencies of polariton modes in an hBN 
spheroid as functions of the aspect ratio A — a±/a z . The 
red lines are exact solutions of eq (SI). The blue circles are 
from the EBK quantization method. The labels are the mode 
indices (Z, m , n) with n — 0 and n > 0 being surface and bulk 
modes, respectively. For the (9,2,0) branch which is classi¬ 
fied as a surface mode, the left part of the blue circles is from 
the surface EBK quantization method, and the right part is 
from the bulk one. The dark green line is from the uniform 
approximation method (Supporting information). 



Aspect Ratio 

FIG. 4. Frequencies of representative periodic orbits as func¬ 
tions of the spheroid aspect ratio. The thick black lines 
are for the bulk orbits B i and B 2 with the period ratios 
t^T 1 : Tq 1 : rf 1 = —2:1:0 and —4:1:0. The thin 
blue lines are for the surface orbits Si, S2, S3 with the period 
ratios rf 1 : t^ 1 = 2 : 1, 1 : 1, and 1:2. The insets show such 
orbits in the real space. For the surface orbits, they include 
all the orbits of the given type passing through the equatorial 
point facing the viewer: two for each Si and S2 and one for 
S 3 . Despite similarity to Figure 3, there is no direct relation 
between the dispersion curves of classical periodic orbits and 
those of quantized eigenmodes, either surface or bulk ones. 


and eq (S16) cannot be satisfied. [In fact, eq (S16) fails 
to give a solution already shortly before £ drops below 
£ c .] The structure of the HSM in this case can be un¬ 
derstood from the following physical picture. The HSM 
must exponentially decrease into the interior of the gran¬ 
ule. It can be viewed as a wave with an imaginary p% 
outgoing from the surface into the bulk, i.e., a surface- 
reflected wave generated in the absence of an incident 
one. Therefore, the amplitude ratio of the two waves is 
formally infinite. On the other hand, this ratio equals to 
e 2<5 , and so the HSM with imaginary p% exists if e lS 00 
or tan | = —i. Using this condition, eq (S59), and ex¬ 
pression for the phase shift S (Supporting information), 
one can numerically solve for the eigenfrequency of the 
HSM for any given l and m. The results of such cal¬ 
culations are illustrated in the left part of the (9,2,0) 
curve in Figure 3. They demonstrate a good agreement 
with the exact dispersion curve at small A where this 
approach is justified. We note that the agreement can be 
greatly improved if the EBK formalism is replaced by the 
so-called uniform approximation, which also enables one 
to smoothly connect small and large A parts of the dis¬ 
persion curve (Supporting information). Lastly, one can 
check that the quantum numbers A = l — m and p = mn 
of the HSM still obey the EBK rules, 

j) p e d0 = -7r(2A - 1), P(f> = p, (15) 

applied now to the effective surface Hamiltonian 


H S =Pg + 



1 


sin 2 £ c 


pI 


(16) 


at energy Hs = 0 . 

While the assumption of a spheroidal granule simpli¬ 
fies the theoretical analysis, one may ask if is it possi¬ 
ble to make some correspondence between such a theory 
and the available experiments that were all done with 
HM samples of non-spheroidal shapes. Our tentative 
answer is as follows. The modes observed in truncated 
hBN nanocones , 3 which were previously called ‘volume- 
confined’ are, in fact, similar to a subset of our HSM, 
specifically, (Z,m, n) = (Z, 0 , 0 ) and (Z, 1,0) modes. The 
modes of cuboidal hyperbolic metamaterials 18 and hBN 
slabs 2,13 are conceptually similar to our bulk modes. 
However, indexing them with Z, ra, orn would be tenuous 
as the conserved quantities in such systems are consider¬ 
ably different from those of prolate spheroids. (For exam¬ 
ple, translational momenta in slabs vs. angular momenta 
in spheroids.) 

Periodic orbits. Classical dynamics can prominently 
impact the structure of quantum energies and quantum 
wavefunctions 25 . In particular, the latter may contain 
‘scars’ — enhanced intensity lines — along these classical 
trajectories. An orbit on an invariant torus 39 defined by a 
set of coordinates i is periodic (closed) if the ratios of the 
individual periods of motion t* are rational numbers. For 
our bulk orbits, the condition is : tq : = z\ : Z 2 : £3 

and for the surface periodic orbits, it is tq : = z\ : 

where all zf s are integers. Figure 4 shows the eigenfre¬ 
quencies of two bulk and three surface periodic orbits as 
functions of the spheroid aspect ratio. We expect that at 
such frequencies the field distribution created by polari- 
tons excited by external sources should exhibit regular 













5 



1400 1450 1500 1550 1600 

Frequency (cm -1 ) 

FIG. 5. Purcell’s factor of a dipole emitter near a hBN 
spheroid of aspect ratio A — tanh 1 = 0.761 and long semi¬ 
axis a z = 500 nm. The emitter is positioned on the x-axis 
(the small yellow arrow). The black and blue lines are com¬ 
puted assuming hBN phonon damping rate T = 4.0cm _1 and 
T = 7.0 cm -1 , respectively. The red line is for T = 7.0 cm -1 
with radiative damping included (see Supporting informa¬ 
tion). The inset depicts the distribution of the absolute value 
of the E x electric field component in the x-z plane at the 
(2,1, 0) resonance frequency. 

geometrical patterns. Below we verify this prediction by 
direct numerical calculations. 

Response to a dipole. A peculiar property of the 
spheroid is that the dipole moment of all m > 1 modes 
exactly vanishes, and so they have extremely weak cou¬ 
pling to far-field radiation. Furthermore, while the dipole 
moment of the bulk, i.e., n > 0 modes is nonzero, it is 
quite small. Detection of all such modes in conventional 
optics experiments 3 will be challenging. However, obser¬ 
vation of these modes may be possible using scanning 
near-field optical microscopy. The latter technique uti¬ 
lizes a sharp metalized tip to perturb and measure the 
system response locally. Crudely, one can model the tip 
as a point dipole and the measured signal as the elec¬ 
tric field created by the system at the location of such a 
dipole, see 40-42 and references therein. The same quan¬ 
tity determines Purcell’s factor — the enhancement of 
the radiative decay of a dipole emitter 43 . 

We assume that the emitter and its dipole moment d 
are in the x-z plane. (Here and below the common fac¬ 
tor e~ luJt is suppressed.) We expand the inner and outer 
potentials in terms of spheroidal harmonics, i.e., the ex¬ 
pressions appearing on the right-hand side of eqs (7a) 
and (7b), cf. Supporting information for details. In an 
ideal lossless HM, the electric field outside (reflected by 
the surface of the spheroid) would diverge at each eigen- 
frequency. If the measured optical constants 3 of hBN are 
used, this divergence is replaced by a finite-width reso¬ 
nance. Purcell’s factor, which is proportional to the field 
induced by the reflected wave at the dipole position, ex¬ 
hibits resonances as well, see Figure 5. The strength 
of the resonances and the frequency spacing between 


them decrease as the indices Z, m, and n increase. As 
a result, low-order resonances give distinct sharp peaks 
while high-order resonances merge into a smoothly vary¬ 
ing background. The latter is similar to the broadband 
Purcell effect near the surface of an infinite HM 44-48 . The 
major resonances are due to the (Z, 1 , 0) modes. Note that 
the perturbing dipole is assumed to have the same am¬ 
plitude d at all uj. However, in the scanning near-field 
experiments instead of such a fixed dipole, one has a po¬ 
larizable tip. The back reaction of the nanogranule on the 
tip is expected to cause a small but observable red shift of 
the resonances. This shift can be modeled using recently 
developed analytical and numerical approaches 40-42 and 
studied experimentally by comparing the far-field spec¬ 
trum of a sparse array of identical granules 3 with the 
near-field spectrum of a single granule. Both types of 
efforts can be subjects of a future work. 

The electric field distribution at sharp resonances is 
dominated by the resonance mode. An example is shown 
in the inset of Figure 5 for (Z,m,n) = (2,1,0). This 
field distribution has the nodal structure of the spherical 
harmonic but shows no ‘scars’. However, ray-like pat¬ 
terns do appear at the periodic orbit frequencies. Fig¬ 
ure S4 B\ depicts the field distribution at the frequency 
of the bulk periodic orbit B\ of Figure 4. The shape 
of the high-intensity ray patterns matches the classical 
periodic orbits (magenta lines). The reason why they 
dominate the field distribution can be understood by 
imagining that it is a superposition of fields created by 
wavepackets launched from a finite-size region facing the 
dipole. Wavepackets whose launch points belong to a 
short periodic orbit create a strongly concentrated elec¬ 
tric field. Other wavepackets follow quasiperiodic classi¬ 
cal trajectories that spread all over the spheroid, giving 
an approximately uniform background. Similar behavior 
is found near the frequency of the periodic orbit B 2l see 
Figure S4 B 2 . 

Near-field imaging experiments are expected to be 
most sensitive to the electric field distribution on the sur¬ 
face of the granule. Panels Si-S 3 of Figure S4 show ex¬ 
amples of such distributions projected on the x-z plane. 
They demonstrate directional ray patterns at the fre¬ 
quencies of the surface periodic orbits S 1 -S 3 (Figure 4). 

Conclusions. We investigated basic properties of con¬ 
fined polariton modes in spheroidal nanogranules of polar 
hyperbolic materials. A physically transparent ray-optics 
method for computing eigenfrequencies and wavepacket 
dynamics of the polaritons was presented and its accu¬ 
racy verified by comparison with the exact analytical re¬ 
sults and numerical simulations. We also suggested how 
to probe these polariton modes experimentally using ex¬ 
ternal dipole sources and/or scanned near-field optical 
microscopy. 

There is a number of other directions to explore. 
For example, we restricted our analysis to the hyper¬ 
bolic regime realized at frequencies inside the hBN Rest- 
strahlen bands. The change of the polariton isofrequency 
surfaces from hyperbolic to elliptical at the extremes of 
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FIG. 6. B i, B<i\ False color plot of \E X \ in a meridional cross 
section of an hBN spheroid due to a vertically polarized dipole 
source just above the north pole. The phonon damping rate 
is r = 7.0cm _1 . Outside the spheroid, the dipole’s own field 
is subtracted away, for clarity. The frequencies in B\ , B 2 are 
1555, 1494 cm -1 , which match the frequencies of B\ and B 2 
in Figure 4 for the chosen aspect ratio A — tanh 1 « 0.761. 
The magenta lines are the bulk periodic orbits. S 1 —S 3 : False 
color plot of Ez at the surface of the spheroid projected onto 
the meridional plane. The dipole is just above the surface at 
the center of each plot. The frequencies in S 1 -S 3 are 1557, 
1535, 1488 cm -1 , same as in Figure 4 for the chosen aspect 
ratio A = tanh 0.5 = 0.462. The magenta lines are the surface 
periodic orbits. 


these bands is a topological transition. One may want 
to investigate signatures of this intriguing transition in, 
e.g., Purcell’s factor. 49 

The studied phenomena may also have far-reaching 
technological implications. One can imagine a whole new 
class of polaritonic devices that would include nanoreso- 
nantors, hyperlenses, infrared photon sources, etc. Such 
devices would be deeply sub-diffractional and low-loss be¬ 
cause phonon-polaritons are immune to electronic losses 
that plague conventional metal-based plasmonics. Our 
general approach may be useful for design and optimiza¬ 
tion of these devices. 

We thank A. Relano and M. Berry for useful dis¬ 
cussions and also F. Guinea for comments on the 
manuscript. 

■ AUTHOR INFORMATION 
Corresponding Author: 









7 


Supporting information 


B. Hamiltonian optics 


I. EIGENMODE DISPERSION 

A. Radial index of the modes 

A nanogranule made of a hyperbolic material possesses 
multiple bulk polariton modes corresponding to the same 
‘angular’ indices l and m. Such modes can be indexed 
with the ‘radial’ quantum number n, as described be¬ 
low. Consider the exact eigenmode equation [eq (9) of 
the main text]: 

lnPr(cos?) = In QY 1 (cosh fj). (SI) 
at ; CLT) 

Following Walker 6 , this equation can be written as 

£ ( l W l_ | /; | V _ 2 _ 

ytan 2 ^ xf tan 2 £ — 1 + 

= — cosh fj — — —— In Q7 1 (cosh fj), 
a cosh 77 

where 0 < xi < 1 are the positive roots of the Legendre 
function P™(x) sorted in the ascending order, N = [(l — 
|m|)/ 2 ] is the number of such roots, [z] is the integer 
part of z, k = 2 N — n + \m\ is either 0 or — 1 , and tan£ 
is defined by eq (5) of the main text: 

tan £ = i — . (S3) 

a z ^/£± 

The right-hand side of eq (S 2 ) is a positive finite number, 
while the left-hand side is a sum of poles that occur at 

tan 2 £ — x~[ 2 — 1 . (S4) 



The approximate eigenmodes of our system can also 
be found by combining the Hamiltonian optics approach 
and the Einstein-Brillouin-Keller (EBK) quantization 
rules 24,30 . In this approach the polariton eigenfunctions 
are zero modes of the effective bulk Hamiltonian 


H b = SijPiPj , (S5) 

which describes the region filled by the hyperbolic 
medium. The eigenfrequency uj is contained implicitly 
in the dielectric tensor £ij(uj). It so happens that the 
boundary condition for <f> can be written in terms of a 
single quantity — the reflection phase shift S — defined 
below [Eq. (S24)]. This fact leads to the existence of 
three conserved quantities in the problem, which implies 
that the system is integrable. Two of such conserved 
quantities are obviously the energy (equal to zero) and 
the z-axis angular momentum L z . The third conserved 
quantity L 12 is introduced shortly below. 

We start with analyzing classical dynamics of the sys¬ 
tem. After a canonical transformation to coordinates 
(£, 0, </>), defined by [eq (4) of the main text] 

p = — iby/e± sin £ sin# , z = b^fe~ z cos £ cos 6 , (S6) 

b 2 = e~ x a 2 - , (S7) 


where 0 < £ < £, £ < 0 < 7r — £, the Hamiltonian H B 
becomes [eq (12) of the main text] 


H b 


(d(x,y,z)\ a / d(x,y,z) 

Kd&O,*) 

sin 2 £ — sin 2 9 sin 2 £ sin 2 9 
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3 


(S8) 


In addition, if m ^ 0, there is another pole at tan 2 £ = 0. 
It is easy to see then that eq (S2) may have multiple 
solutions, as stated above. The number of such solutions 
found in a particular frequency range depends on the 
permittivities £_l(cj) and £ z (uj), which enter tan£. In 
hBN the hyperbolic response occurs in two separate mid- 
infrared bands. The upper band, cj™ < uj < cj^ 0 , is a 
type II HM 8,50 , e± < 0, £ z > 0. As frequency uj changes 
from cjJ 0 to co^ 0 , £_l changes from — oc to 0 while e z is 
positive and approximately constant. To find the number 
of the solutions of eq (S2), one just counts the number 
of the poles crossed by tan 2 £ as frequency changes. We 
can index these solutions by an integer n, which is equal 
to zero if the pole is tan 2 £ = 0 and equal to i if the pole 
originates from Xi, eq (S4). One concludes that n runs 
from 1 to N for m = 0 and from 0 to N for m ^ 0. At 
frequencies that belong to the lower band, hBN behaves 
as a type I HM, 8,50 £_l > 0, e z < 0, and similar analysis 
yields that n runs from 1 to l — \m\ — N. Therefore, the 
total number of the solutions in both bands combined is 
equal to / for m = 0 and l — \m\ + 1 for m ^ 0, as stated 
in the main text. 


The existence of the third conserved quantity L 12 be¬ 
comes evident when one goes through the standard pro¬ 
cedure of separation of variables. In our case, where 
H b = 0, the separated expressions for the momenta are 


Pi 

—i-. 

l 12 - 

Li 

(S9) 

\ 

sin 2 £ ’ 

Pe 

= ± \ 

/ r 10 _ 

Li 

(S10) 

1 3j12 

sin 2 9 ’ 

P<t> 

= L Z 



(Sll) 


The position of the caustic is given by 


£ c = arcsin 



(S12) 


If 0 < £ c < £, where £ is given by eq (S3), the momenta 
p^Pe are real in the rectangular region £ c < £ < £, £ < 
0 < 7r — £, see Figure 2b of the main text. Conversely, 
if £ c exceeds £, no classically accessible region inside the 
spheroid exists, see Figure 2d of the main text. 















Note that the velocity of this fictitious motion has com¬ 
ponents v a = dHB/dp a . Since 0 > £, the signs of vq and 
Pq are the same but the signs of v% and p% are opposite. 

Equations (S9)-(S11) specify a hypersurface in the six¬ 
dimensional phase space that has the topology of a three- 
dimensional torus 39 . According to the EBK rules, the 
total phase acquired across any closed loop on this torus 
should be an integer multiple of 27r. The phase must in¬ 
clude a phase shift of — | upon crossing the caustic and 
the reflection phase shift (s) at the boundary. We found 
(cf. Sec. IF) that due to the sign structure of the ve¬ 
locity components, the vertical boundary segment £ = £ 
(region 2 in Fig. 2 of the main text) and the two hor¬ 
izontal ones 0 = £, tt — £ (regions 1 and 3 in Fig. 2 of 
the main text) have opposite phase shifts, respectively, 
S and —S. Therefore, the EBK quantization conditions 
have the form 


£ 


J 1 \Pt\d£ ~ | +<5 = 2to/, 

(S13) 

£c 


P 


2 / \pe\dO — 25 = 2 tt\ , 

(S14) 

£ 


2ttL z = 27T/i . 

(S15) 


Without loss of generality, the quantum numbers (/i, A, v) 
can be taken to be nonnegative integers. From eq (S15) 
we see than L z = p. We now need to express the remain¬ 
ing classical integral of motion L \2 in terms of the quan¬ 
tum numbers. Using eqs (S9), (S 10 ), (S13) and (S14), 
we obtain 


/ \ - ^7 - 9 + 6 ’ 

J V sin £ c sin £ 2 

£c v 

I - 

2 /i / . / —1-, dO — 25 = 2nX , 

7 V sin 2 £ c sin 2 9 
f 


which can be rewritten as 

2 <pZ(€,€c) ~^+S = 2nv, 
4^ (^,t)-26 = 2nX. 

Here we defined 



(516) 

(517) 

(518) 

(519) 


(5 20 ) 

(5 21 ) 

(5 22 ) 


Compared to eq (24) of the main text, here we add the 
infinitesimal quantities ‘-MO’ in the definitions of A(£, £ c ) 
and H(£, £ c ). These infinitesimal quantities have no effect 
a t £ < £c but they will be important in Sec. II where we 
consider £ > £ c to describe the surface modes. Combining 
these equations, we obtain the expression for £ c : 


sin £ c = 


+ A + fi + ^ 


(S23) 


The final step of the EBK procedure is to account for the 
boundary conditions, which entail a certain equation for 
S. This equation can be written as (cf. Sec. IF) 


S 1 — e lS 

tan - = i -7Y 

2 l + e zS 

/ i i \ 1/2 

i 

\ sin^ £ c sin z £ 


which is the same as eq (16) of the main text. Note 
that fj is determined by the aspect ratio of the spheroid 
a±/a z = tanh fj while £ c is fixed by the quantum numbers 
via eq (S59). For each given set of these parameters, 
the system of equations (S18) and (S24) can be solved 
numerically for cj, the implicit argument of an d e z , to 
find the desired eigenfrequency. 


sin^ £ c sinh 2 77 


1 


1 


(S24) 


C. Correspondence between the EBK and the 
exact eigenmodes 


We are to compare the following two eigenmode equa¬ 
tions. The first one is from the exact solution, eq (SI). 
The second one is from the EBK method, eqs (S18), and 
(S24). To do the comparison, we use the asymptotic 
forms of the associated Legendre functions: 
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P^(cos 0) ~ < 


(cos 2 0 C — cos 2 O) 1 ^ 


cos (<£>f*(0, e c ) - 0 , 0 C < <9 < |, 


1 c m / \ 1 

2 7 2 , ~ 2fl xi/4 exp(-^p(6>,6> e )), — r « 9 < 9 C 

z (COS J 0 — COS^ 0 C ) ' v ' L + 2 


<f = (-i r 


1 (Z + to)! 


_7T l + \ (l - m)\ 


1/2 


L 0 r = 


l + 2 


(S25) 


(S26) 


These expressions are valid for large l and to. They can be derived 51 applying the semiclassical approximation to the 
Legendre differential equation. Similarly, for the Legendre function of the second kind one obtains 


7r r rn 

QT( coshr,)^-----*- 

z (cosh T] — cos 2 0 C 


17 iexp[-x, m (r?, 6 » c )] , 


*TM C ) = 


m 

sin 


cosh 


/ cosh T] \ 

\ COS 6> c / 


— rn sinh 


/ coth?7\ 
\ cot 0 C J 


(527) 

(528) 


The leading contribution to the logarithmic derivatives of the P f 1 and Qf 1 comes from the cosine and the exponential 
terms, respectively. Keeping only these terms, we obtain 


l + l) - 


AlnPrtcosfl.y 

4= In QTicoshff) ~ = -\l ( l + t, ) + 


sin 2 £ 


tan (vr&Oc) -• 0 , 


dp 


dr] 




2 = ' 


sinh 77 


(529) 

(530) 


Substituting these expressions into eq (SI), we see that 
it can be matched with eq (S18) if we set 0 C = £ c , i.e., if 
we make the following correspondence between the two 
sets of integers: 


/ = 2 z/+ A +/x, m — ]i. (S31) 

The numerical demonstration of this agreement is shown 
in Fig. 2 of the main article. 


D. Hyperbolic surface modes 


The momenta inside and outside of the spheroid satisfy 
the equations 


V\-Ve _ V\ 

sin 2 £ — sin 2 0 sin 2 £ sin 2 0 
V% +Pe + v\ 
cosh 2 77 — cos 2 0 sinh 2 77 sin 2 6 


iy/el\fe~zVi = Pv 


(inside), (S32) 
(outside), 

(533) 

(surface). 

(534) 


Hyperbolic surface modes (HSM) correspond to imagi¬ 
nary p£ and p v . Eliminating these variables from the 


equations, we obtain 


Hs = P 2 e 


1 


1 


sin 2 6 sin 2 £ c 


= 0, 


(S35) 


same as eq (17) of the main text. Function Hs(po,P(f),0) 
can be considered an effective surface Hamiltonian for 
polaritons. Here we defined 


. 2 /- ~ -7 \£±£z . 27 ^— 2 - J • (S 36 ) 

sm 2 £ c £±£ z ~ 1 V sm 2 £ smh 2 77 J 

This formula, valid for the surface waves, replaces 
eq (SI 2 ) for the bulk waves. 


E. Uniform approximation 

The lowest-order semiclassical approximation eq (S25) 
diverges near the caustic 0 = 6 C whereas the actual Leg¬ 
endre function P[ n (cos£) remains finite. This is not a 
serious problem for the bulk modes; however, for surface 
waves there is a range of parameters where the bound¬ 
ary £ is close to the caustic £ c . This is the cause of the 
discrepancy between the exact and EBK results seen in 
Fig. 3 of the main text near the aspect ratio A = 1. This 
discrepancy can be greatly reduced by using the uniform 
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2 waves. The gradients of the phases Sjk are the wave 
momenta pj k = — id a Sjk . We assume that only a = z 
and x components are nonzero. The boundary conditions 
are: 

$1 = $2 , £i z d z $i = S2 Z d z <&2 • (S42) 

To satisfy them, all p x k s must be equal. In addition, we 
must have 

All + A 12 ~ A-21 + A 2 2 5 

A-nd z Sn + = — (^ 21^621 + A 22 d z S 2 2 ) • 

Taking advantage of the fact that d z Sn = —d z S 12 and 
d z S 21 = -d z S 2 2 , we get 


Figure SI. Geometry of an auxiliary problem of wave reflec¬ 
tion at the boundary of media with permittivity tensors e\ 
and £ 2 . 



where 


l/l + t 1 — A /A 2 l\ 

2 yl — t 1 + tj \A 22 J ’ 


(S43) 


approximation 52 for the Legendre function: 


Pr(cos< 9 ) ~ V™? 


c 


= 3 


1 V cos 2 0 — cos 2 0 C 
12/3 


1/4 


C = e 


~^no,e c ) 


Ai(C) i (S37) 
(S38) 


£2z dz&n 

£\z d z Su 


(S44) 


Setting A 21 =0, which means there is only outgoing wave 
in medium 2, we get the reflection coefficient 


1 + t 

1 ~t ' 


(S45) 


Here Ai(£) is the Airy function. If we apply this 
approximation to the left hand-side of eq (SI), use 
In Qf 1 (cosh rf) —>q m (? 7 , 0 c ) on the right-hand side, and 
keep the leading terms only, then eq (S34) gets replaced 

by 


i \fzx\fz~z 


1 d( Ai ; (Q 
i d£ Ai(C) 


= l 




dr\ 


r]=r] 


(S39) 


Figure 3 of the main text shows an example of applying 
eq (S39) to computing the (9,2,0) surface mode of an 
hBN spheroid. It yields an excellent agreement with the 
exact eigenfrequency of this mode. 


If the wave in medium 2 is evanescent, i.e., if momenta 
P 2 k are imaginary, then the reflection phase shift S is real. 

Next, we turn to our original problem of internal po- 
lariton reflection at the surface of a suspended nanogran¬ 
ule. The problem can be reduced to the one solved 
above using the special choice of coordinates: coordinates 
(£, 0, 0) inside the granule [eq (S 6 )] and the spheroidal co¬ 
ordinates ( 77 , 0 , </>) in vacuum outside [eq ( 2 ) of the main 
text]. Equation (S42) becomes: 

$1 = $2, iy/ely/e^d^! = dr,$ 2 ■ (S 46 ) 

Therefore, eq. (S45) holds after the following trivial 
change is made: 


F. The phase shift of internal reflections 


1 drqS 2 l _ 1 Pr] 

d ? Sn iy/ely/e^ ' 


(S47) 


To compute the internal reflection coefficient of polari- 
tons off the spheroid surface, we first consider an auxil¬ 
iary problem of reflection at the interface of two media, 
1 and 2, with diagonal dielectric tensors and e 2 , re¬ 
spectively, see Figure SI. To solve this latter problem we 
write the scalar potentials of incident and reflected waves 
as follows: 


$1 = A n e Sl1 + A 12 e Sl2 , (S40) 

$ 2 = A 21 e S21 + A 22 e S22 . (S41) 


Taking advantage of eq (S9) for p% and the similar ex¬ 
pression for momentum p v outside the spheroid, 


Prj = i 



we recover eq (S24). 


(S48) 


II. PERIODIC ORBITS 


The first index in the phases Sjk labels the medium, and To study classical periodic orbits of the polaritons, 
the second distinguishes incident k = 1 and reflected k = it is convenient to perform a canonical transformation 
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Figure S2. (Color online) Left: periodic orbits for different values of ( t £, tq, t </>). The red, green, and blue correspond to, 
respectively, \tqIt<$>\ — 1, 1/3, and 3. Parameter |r^| decreases from left to right. Smaller r {s yield smoother orbits. The value 
of It^/t^I can be inferred from the number of radial bounces that occur during one azimuthal (0) cycle. Center: the dispersions 
(electromagnetic held frequencies) of the periodic orbits as a function of the aspect ratio A = a±/a z of the granule. The same 
color code as in the left panel is used. As |r^| decreases, the distance between adjacent curves of the same color decreases. 
The t t —y 0 limit corresponds to whispering gallery trajectories grazing along the surface of the granule. The first 40 inverse 
integer values of |r^| are included. Right: the periodic orbits dispersions for A = 0.5. The orbits include \tqIt<$>\ — i/j, where 
{i,j} 6 {1, 2, 3}. The values of |r^| are the same as in the central panel. 


to action-angle variables. 39 The new momenta are the 
actions of the three independent loops on the three- 
dimensional torus specified by the constants of motion, 
cf. eqs (S9)-(S11): 



Jo 


J(f) 


p e d0 = 2 





<j> p^dcj) = 2nL z . 


(S49) 


(550) 

(551) 


The motion of each action-angle coordinate pair is peri¬ 
odic with the period 


n = 





(S52) 


As explained in the main text, polariton wavepackets fol¬ 
low the same trajectories in the real space and thus also 
in the space of angle variables as the fictitious particle 
with Hamiltonian Hb • However, the rate of change of 
the angle variables for polariton wavepackets is differ¬ 
ent from dH/dJi by a certain overall factor. Therefore, 
the periods t* themselves do not have a direct physical 
meaning but their ratios do. The phase-space trajectory 
or “orbit” is closed if these ratios are rational numbers. 


Each Ti can be represented by a certain Jacobian. For 
example, is given by 


tT 1 = 


0H b = djH^Je, J+) 

<9J ? d(Jz,Jo,J<j>) 
d{H B ,J6,Jt) L\2, L z ) 


d(H B ,Li2,L z ) d(J^,Je,J^) 


_ d(H B , Jg, J^) 
~ d(H B ,L 12 ,L z ) 

dJe 


-l 


<9(Jg, Jjh J<j>) 

d(H B ,L 12 ,L z ) 

dpi^Je) I" 1 

dLuJ L z ,h b [d(H B ,L 12 )_ • 

For the other two periods, tq and r^, we obtain 

i -i 


-i 


-i 


dJ(: 

~dL~ 


12 


1 Jo) 

2tt d(L 12 ,L z ) 


d(Jo,Jj) 

l z ,h b [d(H B ,L 12 )\ 

1 “I 


d(H B ,L 12 )_ 


(553) 

(554) 

(555) 


Their ratios can be reduced to the following form: 


-l 


: t, 


-l 


• T~ l ~ 

• 1 sh - 


dJ e 

dL 


dJj . 1 d(Jg, Jo) 


12 


8L 12 ' 2tt d(Li2,L z ) ’ 

(S56) 


where the first two derivatives are to be taken at fixed 
L z and fixed Hb = 0. After some algebra, we obtain the 
explicit formulas 

T T 1 : Tp : Tp = 7T - 2A : -A : (B - A) sgn L z , (S57) 
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where A = A(£, £ c ) and B = F>(£, £ c ) are defined by 
eqs. (S21) and (S22). To get a particular periodic or¬ 
bit, we follow these steps. First, we choose the period 
ratios to be desired rational numbers. Next, we deter¬ 
mine A and B consistent with this choice. Next, we solve 
for the constants of motion L z and L 12 from eqs. (S21) 
and (S22). Finally, the orbit is generated and plotted 
using eqs. (S9)-(S11). In general, the orbits can have 
very complicated shapes, as illustrated in Figure S2 (left). 
Roughly speaking, the ratio \t^/tq \ determines the topol¬ 
ogy or the winding number of the orbit whereas \r^\ de¬ 
termines the typical radial distance of the orbit from the 
center of the spheroid. As \t%\ decreases, the orbit is 
pushed closer towards the surface of the spheroid. In 
the limit —»• 0, the orbit becomes a smooth trajectory 
grazing along this surface. This kind of trajectories are 
similar to the whispering gallery modes well known in 
ray optics and acoustics. Therefore, they can be con¬ 
sidered a generalization of the whispering gallery modes 
to the present case of the indefinite Hamiltonian Hb • 
For positive-definite Hamiltonians it has been rigorously 
proven 31,53 that the motion along trajectories sufficiently 
close to the surface of a smooth billiard is regular. There¬ 
fore, such whispering gallery modes are subject to the 
EBK quantization rules. 24 We expect that the same prop¬ 
erty holds for indefinite Hamiltonians as well. 

In fact, a precise relation between classical periodic 
orbits and quantization should exist. According to the 
trace formulas given by Gutzwiller 25 for chaotic Hamil¬ 
tonian systems and by Berry and Tabor 54 for integrable 
ones, the density of states (DOS) of the quantized eigen- 
modes can be represented by a sum over the periodic or¬ 
bits. However, in the present case of the indefinite Hamil- 
toninan, although the ray dynamics in a spheroidal par¬ 
ticle is of course integrable, the density of states (DOS) 
is divergent without a momentum cutoff. Therefore, if 
one carries out the summation in the Berry-Tabor for¬ 
mula 25,54 , one should get infinity not only at some dis¬ 
crete frequencies that are equal to resonance frequencies 
but in fact at all frequencies in the Restshrahlen band. 
How to treat these divergencies is an intriguing problem 
for future work. 

The role of short periodic orbits in our system is also 
unconventional. A fair approximation to the exact DOS 
of a billiard system with the usual quadratic Hamilto¬ 
nians H = p 2 /2m, can be obtained including only the 
contributions of the shortest orbits. 54 However, in our 
case no obvious features of the eigenmode spectra near 
the electromagnetic frequencies uj corresponding to short 
periodic orbits can be identified. We speculate that the 
geometric length of the orbit may not be a relevant quan¬ 
tity for systems with indefinite Hamiltonians such as Hb • 
Note that the frequencies uo of families of orbits hav¬ 
ing the same \t^/tq | and decreasing |r^| tend to clus¬ 
ter together. The corresponding frequencies converge 
to certain value that is a function of the aspect ratio 
A = a_\_/a z , see the central panel in Figure S2. Such uj 
are plotted in Figure S3. Naively, a high density of peri- 
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Figure S3. (Color online). Frequencies of the whispering 
gallery periodic orbits as functions of the aspect ratio A. The 
orbits have period ratios {re/r^l — i/j with {i,j} — {1, 2, 3} 
(same as in Figure S2, right). The inset illustrates the shape 
of the orbits in the real space. 


odic orbits near these whispering galley frequencies may 
lead to enhancement of the properly regularized DOS. 
This problem remains to be understood. 

As discussed in Sec. ID, the momentum p% can be 
imaginary for large enough angular momenta L z or at 
small enough aspect ratios A. Such waves are not the 
whispering gallery waves. Instead, they are HSM de¬ 
scribed by the surface Hamiltonian Hs • The analysis of 
the periodic orbits of the HSM is simpler because there 
is only one period ratio, 


rF 1 : tT 1 = 1 - 


sin £ c 


(S58) 


The HSM orbit is closed if r^ 1 : r^ 1 — n\ : n 2 where n\ 
and 77-2 are integers. Dispersion of several such orbits as 
a function of A are shown in Figure 4 of the main text. 
Comparing eq (S58) with the EBK condition 


1 _ 2A +1 

sin£ c 2/i 


(S59) 


which follows from eq (19) of the main text, we see that 
roughly a quarter of all possible HSM periodic orbits 
(those with odd n\ and even 712 ) are simultaneously EBK 
eigenmodes. 


III. RESPONSE TO A DIPOLE 

A. Quasi-static approximation 

In this section we outline the steps needed to calcu¬ 
late the field created by the nanogranule in response to 
a nearby oscillating electric dipole. We assume that the 
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dipole is located at a point R in the x-z plane. Let the 
spheroidal coordinates of R be (770, # 0 , 0 ) with 770 > V- 
The local direction of the coordinate lines is specified by 
the vectors 


1 cosh 77 sin 9 p + sinh 77 cos 0 z 
a cosh 2 77 — cos 2 6 

1 sinh 77 cos Op — cosh 77 sin 0 z 
cl cosh 2 77 — cos 2 6 


(5 60 ) 

(5 61 ) 


4 >=- 




a = 


A _L > 


(S 62 ) 


a sinh 77 sin 6 ’ 

where p, z, and 0 are the unit vectors in the radial, z, 
and the azimuthal directions, respectively. Suppose the 
dipole moment d is also in the x-z plane, then it can be 
defined in terms of two coefficients, d^ and d#, such that 


d v = (r) d) , d e = (0d) , (S 63 ) 


where fj and 6 are to be evaluated at 77 = 770 and 0 = 6 q. 

As in the main text, we denote by <Fi(r) and ^2(1*) the 
scalar potentials inside and outside the spheroid, respec¬ 
tively. We denote by <Fd(r) the potential of the dipole 
alone. These three potentials admit the expansions in 
series of spheroidal harmonics F z m (0,0): 


Figure S4. (a)-(f): False color plot of \E X \ in a meridional 
cross section of a hBN spheroid due to a dipole source located 
just above the north pole and pointed north. Damping loss 
is neglected. Outside the spheroid, the dipole’s own field is 
subtracted away, for clarity. The frequencies in (a)-(c) are 
1551, 1555, and 1561cm- 1 . In (d)-(f), they are 1490, 1494, 
and 1497 cm -1 . The middle numbers in these sets match the 
frequencies of B\ and B 2 in Figure 4 of the main text for the 
chosen aspect ratio A = tanh 1 0.761. The tilted magenta 

lines run parallel to the polariton group velocity. Si—S 3 : False 
color plot of at the surface of the spheroid projected onto 
the meridional plane. The dipole is just above the surface at 
the center of each plot. The frequencies in Si-S 3 are 1557, 
1535, 1488 cm -1 , same as in Figure 4 of the main text for the 
chosen aspect ratio A = tanh 0.5 ~ 0.462. The tilted magenta 
lines run parallel to the HSM group velocity at the center of 
the image. 


YT = P f (cos 0) COS 7770 , (S64) 

00 l 

*1 = E E (cos orr , $ < ?, (S65> 

l = 0m = 0 
00 l 

= E E D ™ (cosh^ir + , 7? > fj , 

l = 0m = 0 

(S66) 

00 1 

= E E D?P? (cosher , V < Vo- (S67) 

l = 0m = 0 

The expansion coefficients D™ of can be derived from 
the known expansion coefficients 55 of the potential 
of a point charge: 


DT = dd n C lrn = (d v d V0 + d e de 0 ) C T , 


Q m (R) = ^i m (2/ + l) 


(1 — m)\ 


QT( cosh7?o) P™(cos6>o), 


(l + m)\ 

where e m is the Neumann factor: eq = 1 , e m = 2 (m = 1 , 2 , 3 ,...). We obtain 

(l — m)V 2 


(568) 

(5 69 ) 


D m = i m ( 2 l + !) 


_(l + m)!_ 


[dtidrjoQf 1 (cosh Vo)PT ( cos ^o) + deQT ( cos h Vo)de 0 P™( cos ^o)] • (S 70 ) 
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Imposing the boundary condition (S46), we get the matrix equation for the series coefficients: 


( PP(cosO _ -QJ^coshj?) \ = _ f PJ"(cosh rj) \ 

[i^ely/e-^Pr(cos0 -dfiQ?(cosh fj)J [rpJ ~ {r? 1 ) ^P^cosh fj)J 

which has the solution 


j_m 

l l 


^4— [-d fj Q™( cosh??) PP(cosh??) +QP(cosh 7 ?)^PJ n (cosh??)], 

^ om 

detM 1 ’ 


where 


detM = — PJ"(cos£)c^ Q™ (cosh rf) + d^P™ (cos !;) Q™ (cosh 77), 

S™ = — iy/e±^/e^d^P^^os^) PJ n (coshf;) + PJ"(cos^) ^P™(coshf?). 


(571) 

(572) 

(573) 

(574) 

(575) 


When eq (SI) is satisfied, detM vanishes, so that tj n and rj" diverge. This behavior is consistent with having a 
divergent resonant response at the polariton eigenfrequencies. To compute Purcell’s factor (Figure 5 of the main text 
and Figure S5) we need to know the response electric field outside the spheroid. This field is given by 


W 2 = -&[$ 2 (r) - <Mr)] = E^rj + E% e d + , 

oo l 

E 2i = - E E D TrTdi [er(cosh r,)¥{"(&, 4>)}, i*r,,e,4>, 

l = Om = 0 


so that Purcell’s factor is 

f= i+ i^yi im ^ E ^ +deEr ^ • 

Note also that the square of the response electric field is given by 




(E r 2 V ) 


(El 


20 ) 


(Els) 


cosir Tj — cos 2 9 cosn 77 — cos 2 0 


sinh 2 r] sin 2 0 


(576) 

(577) 


(S78) 


(S79) 


The formula for the total field Ei = —9 r 4>i(r) inside the 
spheroid is similar, except it involves coefficients t ™• 

The distribution of the electric field calculated at sev¬ 
eral periodic orbit frequencies of the bulk waves and the 
HSM are shown in Figure 6 of the main text. They 
demonstrate an enhanced amplitude at the locations of 
the classical trajectories launched from the point on a 
surface facing the dipole source. However, at a frequency 
away from the periodic orbit frequencies, wavepackets fol¬ 
low trajectories that spread all over the spheroid, forming 
an irregular background. This effect is most apparent if 
we neglect the damping loss of the media, as shown in 
Figure S4(a)-(f). If the frequency is detuned by roughly 
5 cm -1 , i.e., a mere 0.3% to either side off the periodic 
orbit frequency, the ray patterns disappear. The ray pat¬ 
tern of several surface periodic orbits are also found at 
their frequencies, as shown in panels Si-S 3 of Figure S4. 
If we account for the hBN phonon damping, as we do 
in Figure 6 of the main text, then the polariton prop¬ 
agation length becomes finite. Notably, because of the 
scale-invariance of the problem, this length is not fixed, 


it scales in proportion to the size of the nanogranule. 
This unusual property holds as long as the size of the 
granule is smaller than c/cj, so that the scale-invariant 
quasi-static approximation is valid. Therefore, a better 
measure of damping may not be the propagation length 
but rather the quality factor Q = u/T. In Figure 6 of the 
main text we used the damping rate T = 7 cm -1 , which 
is near the upper end of the experimentally determined 
range . 3 This corresponds to Q ~ 200. Clearly, for such 
T the polaritons still propagate far enough to complete 
the periodic orbits. Additionally, the frequency windows 
for observing these orbits becomes wider, similar to the 
effect of damping on Purcell’s factor resonances. 


B. Radiative correction 

To explain our procedure for computing the radiative 
damping, it is instructive to consider two auxiliary prob¬ 
lems first. We begin with the textbook problem of a 
point-dipole emitter. It is well known 56 that there is a 










15 


small correction to the near field of such a dipole if we 
consider the retardation effect. The correction contains 
a real part, which leads to a shift of the resonance fre¬ 
quency, and an imaginary part, which causes broadening 
of the linewidth and accounts for the energy loss due to 
radiation. We are primarily interested in the radiative 
damping; thus, we retain only the imaginary part: 

E^ = E s taticTi ImE ra d =-y-—— + — &oP- (S 80 ) 

V o 

Next, consider a finite-size nanogranule subject to an ex¬ 
ternal electric field. The field outside is the sum of inci¬ 
dent field and dipole field. They have to satisfy boundary 
condition at the surface of the granule, which leads to 

*op), (S81) 


total dipole moment of the system: the source dipole and 
the spheroid. From the first term of eq (S66), the induced 
dipole moment of the granule has x- and ^-components 

PZ* = Y^ r °i ’ ^nd = 2 ~y d1 A- ( S83 ) 

The potential <f> ra d corresponding to the correction 
i Im E ra d can be written as 

$rad= J2 + Cd')Pl*(coshr ? )y i m , ( S84 ) 

m = 0,1 

where 

c° = —- (—} a 3 , c 1 = —2c° , (S85) 

9 vc/ 

C °S = -J (^) a. c l = ~ C °S ■ (S86) 


2 i 

P = Xo(Eq + i ImE ra d) = xo ( Eq + — 


where xo is the polarization tensor. Solving eq. (S81), we 
get the radiative damping corrected polarization tensor: 


Xo 


i-d(^) 3 xo 


Eq • 


(S82) 


Thus, the right-hand side of eq (S 66 ) changes to 


$2 = + $rad 


oo l 

+ E E d^qt 

l = 0 m = 0 


(S87) 


Finally, let us consider our original problem of a Imposing the boundary condition (S46), we obtain the 
spheroidal granule perturbed by a dipole source. Here the equation for the modified reflection coefficients r z m of the 
radiation correction field should be computed using the dipolar (i.e., m = 0 , 1 ) modes 


M 




d m \ 

WJ 


( P[ n (cosh rj) \ 
ydfj P™ (cosh fj)J ’ 


(S 88 ) 


which can be rewritten as 


ft™) = d+r m —1 ( P P(coshfj) \ 
yF 1 J v s Df 1 ) \dfjPY 1 (cosh. rj)J 

The solution for the reflection coefficient is 


M -c 71 


P™ (cosh rj) \ 
dfjPy 1 (cosh fj J 


s r 


rj = 


det M — c^S] 


i + c: 


DV 


(S89) 


(S90) 


where detM and S™ can be found from eqs (S74) and (S75). If auj/c 1, i.e., if the nanogranule is much smaller 
than the diameter of Wheeler’s radian sphere c/cj, then c m ,c™ 1, and the radiative damping is weak. Far enough 
from the resonances, where r™ is finite, the correction to reflected field to the lowest order in c m ,c™ is 


= Y (cosh r])Y^ + (cosh (S91) 

m = 0,1 

= Y (c m D?r? + Cd^rTQTioosh^Yr + c m D?r?P?{coshn)Y? . (S92) 

rri = 0,1 


The first term in eq (S92) is due to radiative correction 
of the polarizability of the spheroid, the second term is 
the reflected damping field of the source dipole, and the 


third term is the damping field of the induced dipole 
on the spheroid. The corresponding change to Purcell’s 
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Figure S5. Purcell’s factor as a function of frequency. The 
system configuration is the same as in Figure 5 of the main 
text except the long semi-axis of the spheroid is reduced to 
200 nm. The red dots (blue lines) show Purcell’s factor with 
(without) the radiative damping included. 


factor is 


Calculations done according to the above formulas are 
shown in Figure 5 of the main text for the case of spheroid 
with the long semi-axis a z = 500 nm and in Figure S5 for 
a z = 200 nm. The most noticeable effect in these Fig¬ 
ures is the broadening of the resonance peaks. There are 
two sources of such broadening. The first is the intrinsic 
loss of the medium, described by the phonon damping 
rate T. It influences all the modes, i.e., we can see that 
all the resonant peaks become broader as we increase T. 
The second is the radiative damping effect, which broad¬ 
ens the mm 0 and 1 dipolar modes but does not change 
much the linewidths of the remaining (2, 2, 0) and (3, 2, 0) 
modes. However, for a z = 200 nm spheroid, c 1 is already 
so small that the radiative damping correction is negligi¬ 
ble, cf. Fig. S5. 

IV. DIELECTRIC FUNCTION OF hBN 

Although our theory is developed for a general hyper¬ 
bolic material, all the figures are calculated for the polar 
insulator hBN. The dielectric tensor components of hBN 
have the following form: 


Sf = ^i h ( S93 ) £iH=£i(oo) 


K LO ) 2 - K TO ) 2 


K TO ) 2 


■ 


• icjR- 


(S96) 


which is meant to be evaluated at 77 = 770 and 0 = 6$. 
In particular, if the dipole is polarized in the x-direction 
and is located on the positive-z semi-axis, then eq (S93) 
simplifies to 


5f = Re 




(S94) 


where the argument of both Q\ and Pi is equal to cosh 770 - 
Note that Sf is free of the small parameter auj/c. Hence, 
the ratio of the radiative correction and the uncorrected 
Purcell’s factor scales linearly with (acj/c) 3 <C 1. 

At the resonances, rf 1 diverges, so one should avoid 
expanding the denominator in eq (S90). This leads to 
the more accurate formula 


Sf = Re < 


1 - c x r\ 


c 1 + 2 — +T 1 (— 

PI + 1 \P} 


(S95) 


instead of eq (S94). This equation indicates that the 
resonances of the radiation-corrected Purcell’s factor are 
shifted from points r\ — 00 to points r\ = 1/c 1 1. 

Existence of such a shift is expected on general grounds. 
However, we do not attempt to evaluate it because a more 
important contribution to this shift should come from 
the real part of the radiative reaction field neglected in 
eq (S80). 


where i = _L or z and 3 

= 1360 cm -1 , l = 1614 cm -1 , (S97) 

cjJ 0 = 760 cm -1 , co^° = 825 cm -1 , (S98) 

£_l(oc) =4.90, £ 2 ( 00 ) = 2.95 . (S99) 

The results shown in Figure 5 of the main text are calcu¬ 
lated for two values of the damping rate Y = Y z , namely, 
4 cm -1 and 7 cm -1 to illustrate the effect of dielectric 
losses. 

We treated the nanogranule as a continuum medium. 
As observed in recent experiments 2 even for hBN as thin 
as three atomic layers the continuum medium treatment 
that uses the bulk dielectric tensor yields an excellent 
agreement with the observed mode spectra. Validity of 
the continuum medium treatment for few-layer hBN can 
also be justified theoretically based on the phonon dis¬ 
persions of a few-layer hBN calculated by diagonaliza- 
tion of the full dynamical matrix. 57 One can easily check 
that these results match with the continuum medium 
treatment for phonon momenta much smaller than the 
inverse lattice constant, which are relevant for our con¬ 
sideration. The main difference between the continuum- 
medium and microscopic theories is the total number of 
phonon-polariton modes. This number is finite and is 
proportional to the total number of layers. 57 We have 
in mind nanogranules which contain hundreds or even 
thousands of layers. For such granules, our continuum- 
medium theory should be fully valid. 
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